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ABSTRACT 

We report on the first fully consistent conventional cluster simulation which includes terms up 
to post 5 / 2 Newtonian in the potential of the massive body. Numerical problems for treating 
extremely energetic binaries orbiting a single massive object are circumvented by employ- 
ing the special "wheel-spoke" regularization method of Zare (1974) which has not been used 
in large-TV simulations before. Idealized models containing N = 1 x 10 5 particles of mass 
1 Mq with a central black hole of 300 M Q have been studied on GRAPE-type computers. An 
initial half-mass radius of r n ~ 0.1 pc is sufficiently small to yield examples of relativistic 
coalescence. This is achieved by significant binary shrinkage within a density cusp environ- 
ment, followed by the generation of extremely high eccentricities which are induced by Kozai 
(1962) cycles and/or resonant relaxation. More realistic models with white dwarfs and ten 
times larger half-mass radii also show evidence of GR effects before disruption. Experimen- 
tation with the post-Newtonian terms suggests that reducing the time-scales for activating the 
different orders progressively may be justified for obtaining qualitatively correct solutions 
without aiming for precise predictions of the final gravitational radiation wave form. The re- 
sults obtained suggest that the standard loss-cone arguments underestimate the swallowing 
rate in globular clusters containing a central black hole. 
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1 INTRODUCTION 

TV-body simulations of stellar systems containing one or more 
massive bodies have become popular with a number of studies 
(Milosavlevic & Merritt 2001, Baumgardt, Makino & Ebisuzaki 

2004, Baumgardt, Hopman, Portegies Zwart & Makino 2006). This 
type of work has been inspired by recent determinations of massive 
black holes (BHs) which reside in galaxies and intermediate mass 
black holes (IMBHs) have also been considered, particularly in re- 
lation to globular clusters. Although the presence of BHs has mo- 
tivated such work, few papers so far have included the relativistic 
effects that would be required for an appropriate treatment of de- 
viations from Newtonian motions. One reason for this situation is 
the numerical challenge of dealing with point-mass dynamics, both 
as regards the large range in time-scales as well as nearly singular 
interactions involving superhard binaries. Regularization methods 
have traditionally been used to overcome problems of this kind. 
However, existing codes are not tailor-made for dealing with large 
mass ratios and some attempts to employ two-body (Milosavlevic 
& Merritt 2001) and chain regularization (Szell, Merritt & Mikkola 

2005, Amaro-Seoane & Freitag 2006) only succeeded in obtaining 
modest shrinkage of the dominant binary. 

Binary shrinkage by itself is unlikely to be be sufficient in 
achieving relativistic stages for all but the most extreme initial con- 
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ditions. Thus, given two BH binary mass components tubh = 
N 1 ' 2 !^, where rh is the average mass of N field stars, the fi- 
nal semi-major axis for absorbing half the total energy would be 
Ar-^/N in terms of the half-mass radius r\. Such evolution is 
well within practical reach for N = 10 J with present methods. 
However, this would fall well short of the required shrinkage for 
relativistic conditions using conventional length-scales. The ques- 
tion therefore arises whether large eccentricities may be generated, 
thereby reducing the relativistic time-scale significantly by the clas- 
sical factor (1 - e 2 ) 7/2 . 

Results obtained so far have yielded conflicting information 
about the eccentricity evolution but clues are emerging that this 
may well depend on initial conditions. An early simulation based 
on two merging clusters, each containing a massive single object 
inside a cusp-like system of N ~ 1 x 10 5 particles, produced suf- 
ficiently large values (e max > 0.995) for the GR coalescence con- 
dition (Aarseth 2003b). This work employed the time-transformed 
leapfrog method (ttl) which was specially designed to deal with 
a massive binary BH (Mikkola & Aarseth 2002) and included the 
three first post-Newtonian terms. Large eccentricities were also re- 
produced recently by placing two massive objects inside a single 
rotating cluster (Berczik et al. 2006). A new study of three massive 
BHs by direct integration (Iwasawa, Funato & Makino 2006) did 
achieve coalescence by including relativistic energy loss, albeit for 
large mass ratios and short times. 

Although the binary BH problem is rightly receiving con- 
siderable attention, the simpler case of just one massive object is 
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still waiting to be explored fully. So far only one major investiga- 
tion appears to have been directed towards this goal (Preto, Mer- 
ritt & Spurzem 2004). In this study, the attention was focused on 
the growth of a density cusp around a massive central object for 
comparison with Fokker-Planck solutions. In spite of using chain 
regularization (Mikkola & Aarseth 1993), the maximum growth in 
1/dBH was relatively modest, mainly due to the mass ratio of about 
1000 exceeding somewhat the above-mentioned factor of N 1 ^' 2 for 
N in the range 1 - 2.5 x 10 5 . 

On the theoretical side, some efforts have been made to look 
for enhanced relaxation in connection with the loss-cone problem. 
The interesting concept of resonant relaxation (Rauch & Tremaine 
1996) has yet to be tested in large-scale simulations. More recent 
work (Hopman & Alexander 2006) has hinted at a connection with 
the Kozai mechanism. On the basis of Fokker-Planck simulations it 
is suggested that coherent torques change the angular momenta and 
rotate the orbital inclinations, giving an increased tidal disruption 
rate for the central region. 

The present investigation concentrates on the evolution of the 
central subsystem using a powerful new implementation. The ba- 
sic method is a generalization of three-body regularization (Aarseth 
& Zare 1974) which is of the same vintage (Zare 1974) but apart 
from one application to four-body scattering (Alexander 1986) has 
apparently remained dormant. Because of the analogy with spokes 
on a cart-wheel being connected to the hub, it has been termed 
"wheel-spoke regularization" (Aarseth 2003a) since each spoke 
represents mutual interactions with respect to the central hub or 
reference body. Experience with the three-body regularization code 
suggests that the extension to arbitrary memberships would make 
an ideal tool for studying the single BH problem, as was already 
suggested when the phrase was coined. Such a scheme has the 
advantage of the massive object acting as a permanent reference 
body, thereby allowing considerable simplifications compared to 
the standard chain algorithm. At the same time, post-Newtonian 
terms are included as perturbations in anticipation of favourable 
developments for which very high eccentricities would be needed. 
A full post-Newtonian treatment has already been tried in a dif- 
ferent context (Kupi, Amaro-Seoane & Spurzem 2006) for an ex- 
tremely compact system of 1000 point-mass particles with rms ve- 
locity about one percent of the speed of light. The aim here was 
to study run-away evolution by capture and mergers due to gravi- 
tational radiation using two-body regularization, as was first done 
a long time ago with the code nbodys (Lee 1993) without the 
precession terms. 

This paper begins by outlining the basic method and its imple- 
mentation in an A-body system, followed by a description of the 
initial conditions for an equilibrium distribution. We report on the 
results of several similar idealized models which illustrate the ca- 
pabilities of the method. Realistic models which include finite-size 
effects involving white dwarfs are also considered. Some features 
of the post-Newtonian implementation are described, with empha- 
sis on the experimental nature of the scheme. Finally a brief discus- 
sion is given, together with some suggestions for the future. 



2 WHEEL-SPOKE IMPLEMENTATION 

We begin by defining a subsystem of n single particles of mass ra» 
and one dominant body denoted by i = with mass mo. The initial 
coordinates and momenta in the local centre of mass are given by 
q i; pi for i = 0, . . . , n. Let us take p = for the reference body 
and introduce relative coordinates qi = qi q with respect to 



mo. The Hamiltonian then takes the form 

x — > p, 1 v — > T v — > m,i \ — > m,m, 

i—i i<j i=i i<j 

with fii = m.imo/(m; + mo) and Ri — |qi | . 

Canonical variables Qi, Pi with a separable generating func- 
tion are now introduced for each two-body pair im , m by (Zare 
1974) 

n 

w( Pi ,Qi) = J2pT -WQi), (2) 

i=l 

where fi(Qj) connects physical and regularized coordinates for 
spoke index i. The corresponding regularized momenta take the 
form 

P I= A lPl , (i = l,...,n), (3) 

where A; is twice the transpose 4x4 Levi-Civita matrix. Inverse 
transformations yield the relative coordinates and momenta 

qi^A^Qi, p, = i^i, (4) 
l 4 Hi 

from which the final physical variables are readily determined. 
Regularized equations of motion for each spoke interaction are 
finally obtained by introducing a time transformation. In accor- 
dance with standard practice in multiple regularization (Alexander 
1986, Mikkola & Aarseth 1993) we adopt the inverse Lagrangian, 
t' = 1/L, which has proved effective. Differentiation of the regu- 
larized Hamiltonian T* = t'(H — Eq) (where Eo is the total en- 
ergy) with respect to Q and P then yields the equations of motion 
to be integrated. 

The implementation of the wheel-spoke scheme into a gen- 
eral A-body code requires many special procedures. Quite a few of 
these can be taken over from other large simulation codes (Aarseth 
2003a) when due allowance is made for differences in the data 
structure. In the following we summarize some of the most relevant 
features which have much in common with chain regularization. 

To begin with, a suitable subsystem must first be chosen for 
special treatment. The idea here is to select an energetic binary con- 
taining the massive BH surrounded only by a small number of per- 
turbers in the central density cusp. The latter requirement invariably 
implies a relatively small semi-major axis, say obh — R c u where 
R c \ ~ 10r h /N is the adopted close encounter distance for stan- 
dard regularization treatment. Upon selection of a suitable binary, 
a few nearby perturbers are added to the new subsystem which is 
initialized as a composite particle for direct A-body integration in 
the usual way (Aarseth 2003a). 

Subsequently a variety of heuristic conditions have been em- 
ployed to control the subsystem size. Particular attention is di- 
rected towards limiting the maximum membership because of the 
n(n — l)/2 force terms which are accumulated by the accurate 
but time-consuming Bulirsch-Stoer (1966) integrator. Thus the so- 
lution method is based on shared but variable time-steps and each 
step requires many function evaluations for a tolerance of 10~ 12 . 
Likewise, internal members moving outside a specified distance 
are removed from the subsystem and play the role as external per- 
turbers while tidal effects are important. The processes of absorb- 
ing or emitting subsystem members entails careful updating of the 
energy budget in order to maintain conservation at a high level. 

Particles outside the subsystem also play a role in driving 
the evolution. External perturbers are selected for distances d < 
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[2mj/mo7o] 1//3 ^grav, where 70 = 1CP 6 is a small dimension- 
less parameter measuring the relative perturbation at the boundary 
and i? g rav is the gravitational radius (2<jbh for a binary). Since 
the wheel-spoke formulation only yields regular equations for in- 
teractions with the massive body, it is necessary to include a small 
force softening in the singular terms, thereby permitting the smooth 
treatment of near-collisions. The softening length was first taken as 
e = eo-Rgrav with eo — 0.01 and later reduced to eo = 0.001 
without experiencing numerical problems. It is updated every time 
the subsystem is reduced to a dominant binary which occurs fre- 
quently. Consistent potential energy corrections are then made at 
each change in membership. It should be emphasized that soften- 
ing is only applied between the few non-spoke internal interactions 
so that the essential dynamics is maintained. 

The subsystem solution is combined with the standard Her- 
mite block-step scheme (Makino 1991) which ensures synchro- 
nization and consistency. Accordingly, the force on any subsystem 
perturber is evaluated carefully by summation over all the internal 
members. The adaptation to special-purpose GRAPE computers is 
based on the standard nbody4 code (Aarseth 2003a) where chain 
regularization has been bypassed. We note here that possible nu- 
merical problems for large mass ratios in chain regularization are 
circumvented by the present Hamiltonian formulation where the 
different relative motions are treated on an equal footing. Likewise, 
the absence of switching and re-initializing the chain for the same 
membership is beneficial. 

A special feature in the present investigation is the addition of 
post-Newtonian terms in the relativistic regime. This treatment ne- 
cessitates the introduction of physical units. Scaling to total energy 
E = — 1/4 for the gravitational constant of unity and J2 mi = 1 
gives an rms velocity of 2 1 / 2 /2 and velocity unit V* in km s _1 
once the total mass and length scale is specified. Hence the speed 
of light is given by c = 3 x 10° /V in model units. The result- 
ing equation of motion can be written in a convenient form by 
(Blanche! & Iyer 2003, Mora & Will 2003) 

d 2 r _ M 

dt 2 ~ V 2 

where the scaled quantities A and B represent the accumulated rel- 
ativistic effects associated with the separation vector r and relative 
velocity v. The latter is readily obtained from the momentum trans- 
formations when due allowance is made for the non-zero momen- 
tum of mo in the local frame. 

In the present formulation the perturbing force is required for 
consistency with the equations of motion (cf. Mikkola & Aarseth 
1993). Consequently, the desired expression for the dominant two- 
body motion takes the form 

^ + ^ + ^)-+ 

C Z pi r 1 r 

(6) 



[{-l + A) r - + Bv\ 



(5) 



F G r = 



mono 



c 2 c 3 



In view of the increasing number of operations involved for 
the higher orders, it is of interest to implement an efficiency scheme 
depending on time-scale and order. Such an experimental approach 
aims at providing a qualitative description of the inspiralling pro- 
cess without any predictions for the final coalescence event. 

The classical time-scale for radiation energy loss is employed 
for decision-making purposes. In A-body units the instantaneous 
value, a/a, for a large mass ratio is given by (Peters 1964) 



TGR 



%Amira\ g(e) 



(7) 



where g(e) is a known function (~ 4.5 for e = 1). Provided 
the time-scale falls below a specified value, the radiation terms 
A5/2 and B5/2 are added to any Newtonian perturbations for the 
dominant interaction. Other terms are activated progressively on 
subsequent reduction of tgr. Moreover, comparisons with the 
full scheme for a nearly isolated binary with circular or eccen- 
tric orbit yield final coalescence times in essential agreement. En- 
ergy conservation is monitored by separate integration of the post- 
Newtonian perturbation according to 



AE G r = 



v dt . 



(8) 



converted to the appropriate regularized form (Mikkola & Aarseth 
1993). Note that satisfactory conservation does not guarantee the 
solution here and is merely a consistency check; hence comparison 
with known two-body solutions should be carried out. 

Since the present full formulation is applicable for quite large 
deviations from Newtonian dynamics we define GR coalescence in 
the traditional way as three Schwarzschild radii by 

6 (mi + mo) 



^"coal 



(9) 



or about 8 x 10~ in A-body units for the standard parameters. 
This represents over five orders of magnitude in distance range with 
respect to the initial close encounter length scale, R c \ . 



3 INITIAL CONDITIONS 

Each simulation project requires careful considerations concern- 
ing the initial conditions. In the present work, the introduction of a 
massive central body may have a significant effect on other nearby 
particles. Since the emphasis is on studying the evolution of the in- 
nermost part, it is desirable to start with an approximate equilibrium 
distribution. Following an earlier formulation (Aarseth 2003b), we 
adopt a cusp-like stellar density profile 



p(r) = 



r l/2( 1+r 5/2) 



(10) 



Given a central body of mass mo, the corresponding ID velocity 
dispersion is generated by (Zhao 1996) 



a 2 (r) 



P( 



l_ r 



^-[m{r)+m ]dr, (11) 



where m(r) is the enclosed mass within r. 

Several models with N — 1 x 10° equal-mass particles of 
1 M Q are studied. Based on the discussion above, we take the BH 
mass as A 1//2 m which corresponds to mo = 3xl0~ 3 in A-body 
units with rm = 1 x 10~ 5 . This mass ratio is relatively modest 
compared to typical values of tubh — 0.01 used in many investi- 
gations of massive binary components. 

The choice of total mass and half-mass radius determines the 
scaled speed of light in the simulation and hence the influence of 
any post-Newtonian interactions. Since this investigation is a first 
attempt to evaluate the code performance, standard star cluster pa- 
rameters would appear to be outside the likely range of interest un- 
less exceptionally large eccentricities are reached. For this reason 
we first adopt a more conservative length scale of 0.1 pc which 
yields an initial half-mass radius r h = 0.094 pc. With the ve- 
locity scaling V* ~ 66 kins -1 , this makes c ~ 4600 in N- 
body units. Consequently, an energetic binary with semi-major axis 
a ~ 3 x 10~ 5 would need an eccentricity exceeding 0.99 for the ra- 
diation time-scale to fall below 400 A-body units which might be 
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Figure 1. Semi-major axis and eccentricity as functions of time, model G2. 
GR coalescence occurred at t = 92, 104, 117, 120, 130, 144, 152. 
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Figure 2. Semi-major axis and eccentricity as functions of time, model G3. 
Coalescence took place at t = 22, 30, 50, 84, 104, 186. 



a reasonable criterion for initiating the post-Newtonian treatment 
since most simulation times are considerably shorter. This com- 
pares to the smallest initial central binary of size a — 2.5 x 1CP 4 
in a typical model and only one more below 1 x 10 . Hence a sig- 
nificant evolution is required to achieve a factor of 10 shrinkage of 
the most energetic binary, especially bearing in mind the increased 
central velocity dispersion associated with the massive body. 

The early cluster evolution exhibits small departures from 
overall equilibrium as defined by the virial energy ratio. Likewise, 
the number of particles inside the innermost fixed radii do not show 
significant changes, in accordance with expectations for an equilib- 
rium model. 

Several models assume that the interactions are between point- 
mass particles. However, it is also of interest to study finite-size ob- 
jects. The possibility that stars are disrupted by the BH adds another 
complication. We introduce the tidal disruption distance (Magor- 
rian & Tremaine 1999), 



n = 



m 
mi 



1/3 



(12) 



and adopt r* = 5 x 10 _s au for white dwarfs of 1 M . In these 
models, any star inside this distance is removed from the calcula- 
tion and its mass added to the BH instantaneously. The question 
of whether such orbits may be modified by post-Newtonian effects 
before disruption depends on the cluster parameters. 



4 NUMERICAL RESULTS 

We first present some results for idealized models which have been 
studied over at most a few hundred iV-body time units and is suffi- 
cient to illustrate the general behaviour. Figure 1 displays the semi- 
major axis of the innermost bound orbit in model G2 (rh = 0.09 
pc), together with the corresponding eccentricity. Since there are 
seven GR coalescence events, the downwards trend of the semi- 
major axis is replaced by the next most strongly bound orbit, where- 
upon the shrinkage continues. Inspection of the eccentricity graph 
reveals a number of spikes, with both high and low values which are 
the hall-mark of Kozai cycles. The growth in eccentricity preceding 
significant GR radiation loss can be substantial, with pre-relativistic 
values exceeding e ~ 0.99999 on several occasions. Likewise, the 
associated predicted maxima are often very close to unity when the 
inclination is near 90°. However, not all such configurations are 
sufficiently stable for GR conditions to develop. 

The second model (G3), shown in Fig. 2, exhibits a different 
behaviour of the smallest semi-major axis, this time with six in- 
spiralling events. Note that during two epochs (t ~ 50 — 70 and 
86 — 98) the central subsystem is too weakly bound for the spe- 
cial treatment. After an irregular early phase there is a downward 
trend which is barely affected by the last coalescence. At termi- 
nation the semi-major axis decreased below the last plotting point 
to 1.2 x 10~ 6 and would soon have satisfied the coalescence cri- 
terion. During the non-relativistic stage the innermost semi-major 
axis tends to shrink as a consequence of dynamical evolution. How- 
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Figure 3. Semi-major axis and eccentricity as functions of time, model Wl . 
Disruption occurred at t = 25.0, 62.8, 70.2, 88.2, 95.2, 102.5, 104.8. 
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Figure 4. Semi-major axis and eccentricity as functions of time, model W2. 
In this model there was just one disruption at the end. 



ever, the shrinkage would progress much further but for the large 
eccentricity which therefore places a lower limit on the period be- 
fore the GR effect intervenes. Each coalescence is associated with 
a substantial loss of radiation energy with an accumulated value 
AEcn — —181, consistent with the final Newtonian binding en- 
ergies. This may be compared with a small systematic drift of 
—2 x 10~ 5 units in the total energy which should be conserved. 

These early models employ the point-mass assumption as in 
most current work. Based on the promising experience gained for 
point-mass particles it is worth while to enlarge the investigation 
and include finite-size objects. Given the uncertain status of neu- 
tron star populations in globular clusters, it may be more appropri- 
ate to consider white dwarfs which should be present with signif- 
icant abundance in the cores. Since the code can handle a general 
stellar distribution at a specified age, a study of more realistic mass 
functions will be undertaken in due course. 

In the second set of models, the field stars are taken to be 
white dwarfs with the characteristic radii and masses quoted above. 
For continuity we begin with a similar cluster half-mass radius of 
0.09 pc (model Wl). Figure 3 gives the smallest semi-major axis 
as a function of time, together with the corresponding eccentric- 
ity. There are seven disruptions in this interval, with the relevant 
times quoted in the figure caption. All but one of these events show 
significant GR energy loss while the exception is a head-on col- 
lision with large pre-relativistic eccentricity, e = 0.999999. Af- 
ter the seventh disruption, the accumulated relativistic energy was 



AEgh ~ —0.31 which exceeds the initial cluster energy. Once 
again the eccentricity plot reveals the characteristic spikes associ- 
ated with Kozai cycles. Note that the corresponding time-scales are 
sometimes shorter than the plotting intervals so that fine structure 
is often missing. 

The cluster half-mass radius was also varied in order to ex- 
plore the length-scale dependence. First, a large half-mass radius 
with rh = 2.8 pc was chosen for model W2. This gives a smaller 
disruption radius of r t ~ 5 x 10~ 10 compared to 1.6 x 10~ 8 above. 
The results shown in Fig. 4 contain just one disruption at a compa- 
rable time to Fig. 3 (t{ ~ 95). Again a Kozai cycle was apparent, 
with e max ~ 0.9998. A rather long period of circularization then 
took place before termination (e = 0.1, r = 5 x 10~ ). Until this 
time, the previously accumulated relativistic energy loss was small 
which indicates that no other critical approaches were present. 

In the final model W3 a cluster of intermediate half-mass ra- 
dius r h = 0.9 pc was studied. A collision with small predicted 
pericentre for a wide orbit (a ~ 0.02, 1 — e = 10~ 7 ) occurred dur- 
ing the early stage before the compact subsystem formation. All the 
five disruption events exhibited maximum eccentricities exceeding 
0.9999 shortly before termination, evaluated at a non-relativistic 
distance r ~ 10 - 5 . Moreover, in each case the energy loss by 
equation (8) yielded an amount AEgr — —0.5. This shrinkage 
corresponds to a semi-major axis ogr ~ 3 x 10~ 8 which agrees 
with the Newtonian value obtained just before the final plunge. 

The density distribution in the central region is also of in- 
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Figure 5. Semi-major axis and eccentricity as functions of time, model W3. 
Five disruption events occurred at t = 38.7, 75.0, 108.6, 115.6, 143.8. 



terest. At typical advanced stages of inner binary evolution (a ~ 
2 x 10~ 5 ), there are only about 10 members inside a distance of 
300a. Moreover, the second innermost semi-major axis tends to be 
well separated from the innermost one, thereby giving opportuni- 
ties for favourable Kozai cycles. Only a few of the dozen or so most 
extreme eccentricity excursions result in significant energy loss and 
may therefore be considered rare events. 



5 POST-NEWTONIAN ASPECTS 

The results presented above show examples of relativistic effects 
leading to coalescence or disruption. Although the inner binary 
shrinkage is often considerable, the eccentricity plays a key role 
in initiating the relativistic stage. The process by which large ec- 
centricities are attained is mainly due to third-body secular pertur- 
bations known as Kozai (1962) cycles. This requires the perturber 
to remain in a long-lived orbit around the innermost binary with 
initial inclination exceeding about 40°. Conservation of angular 
momentum yields a relation connecting the inclination and inner 
eccentricity in the form cos 2 8(l — e 2 ) = const. The period for the 
induced inner eccentricity variation is given by 



Txozai — 



T ou t ( 1 + '/out \ r, 2 \3/2 r I i\ / n \ 

-7TT- (1 - e ut) ' /(ein, Win, V0 , (13) 

-tin \ 9out / 



where g ut = milirnx + mo) is the mass ratio for the outer or- 
bit with elements a out , e out and the function / which depends on 
orbital elements is usually of order unity (Heggie 1996). 

An expression is also available for the maximum eccentricity 
of the inner orbit, e max , which is of considerable interest (Heggie 
1996). Since the period ratio Tout/Tin may be taken of order 10 
for a long-lived system here, this gives Tkozai — 3000T ou t(l — 
e out) 3 ^ 2 - Hence a Kozai period of about 1000 T out may suffice for 
e out ~ 0.7 and a moderately high inclination, which is frequently 
seen in such simulations. 

A new decision-making scheme has been developed to intro- 
duce the different post-Newtonian terms for separate time-scales. 
Although these terms are increasing in complexity, this is done only 
partly for computational efficiency. The idea of considering the ad- 
ditional perturbations progressively (Aarseth 2003b) has also been 
tested by a stand-alone code for the three-body problerrQ (Aarseth 
& Zare 1974). In the present work, the emphasis is on achieving 
a qualitatively correct description of the evolution towards coales- 
cence or tidal disruption without aiming for detailed predictions. 

Since the main decision-making is based on the two-body el- 
ements a and e, it is desirable to use relativistic definitions when 
appropriate. This procedure was not adopted for most of the work 
reported here. However, the elements are usually evaluated outside 
the semi-major axis of predominantly highly eccentric orbits where 
GR effects are negligible. For a given level of post-Newtonian im- 
plementation, we add the relevant contributions from the two first 
A and B terms of equation (6) according to the relativistic expan- 
sion (Mora & Will 2003). Without these corrections, the classical 
elements show a range of values depending on orbital phase and 
v/c. Now the modified two-body energy and angular momentum 
yield consistent values of the osculating semi-major axis and ec- 
centricity, as can be verified by the available three-body code for 
an isolated binary. We remark that the relativistic contributions 
for the white dwarf models are fairly modest, with typical values 
»/c~4x 10~ 4 for r > a near disruption in the case r n = 0.9 pc. 

In the case of a binary we define appropriate time-scales for 
activating the perturbations in terms of the instantaneous radiation 
time-scale, tgh .. Thus the radiation term 2.5PN is already included 
for tgh, < 1000 TV-body units. The subsequent terms 1PN, 2PN 
and 3PN are then activated below the experimental values of 100, 
50 and 10, respectively. Likewise, the relevant terms are excluded 
for increasing time-scales up to tgr > 2000 for the GR radia- 
tion, with slightly more generous limits being applied in order to 
prevent repeated switching. Three-body tests show that the final 
coalescence times for close binaries of different eccentricity are 
in essential agreement with the full scheme when criteria of this 
type are employed for a range of parameters. Regarding the choice 
of boundary values for different levels, the nominal time-scale ex- 
ceeds the actual coalescence time for circular orbits by a small fac- 
tor depending on c due to the loss of angular momentum. At some 
advanced stage, it may be justified to adopt the unperturbed two- 
body approximation (Peters 1964) or even define premature termi- 
nation in order to speed up the calculation. 

The possible presence of Kozai cycles adds a further compli- 
cation since the precession terms act to de-tune the resonance. In 
particular, the classical "Mercury" precession gives rise to a peri- 
centre advance in radians per orbit, 



See the public three-body regularization code TRIPLE3 at 
|http : / /www . ast . cam . ac . uk /^sverre /mult ireg| . 
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Figure 6. Eccentricity as function of time for three values of c. The elements 
of the hierarchical three-body configuration are defined in the text. 
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which may be substantial for some conditions of interest. Unless 
existing post-Newtonian perturbations are present due to a suffi- 
ciently short radiation time-scale, the first and second-order pre- 
cessions are introduced for Tkozai < 10. Thus for an eccentric 
binary with e = 0.99 the precession condition gives Aui ~ 0.01 
for a typical inner binary of size a = 1 x 10 _J which is usu- 
ally covered by either of these criteria. More extreme stages may 
give rise to a faster advance which is reflected in a shortening of 
the time-steps. The basic three-body code, which contains identi- 
cal post-Newtonian expressions to the full code, was used to com- 
pare the pericentre advance in the case of strong relativistic effects 
while suppressing the higher order terms. Here we employed an al- 
gorithm measuring the angle between successive apocentre turning 
points which gave good agreement with equation (14) and therefore 
provides an independent check on the first-order solution. 

Two small bodies orbiting a much heavier one in relative iso- 
lation may be sufficiently stable for Kozai cycles to be induced. 
All that is required is that an outer orbit increases its energy and 
has a favourable inclination for large eccentricity growth. As has 
been noted above, this process is seen at later stages when the in- 
ner core has been partially depleted. Hence it appears that the fi- 
nal orbital shrinkage before relativistic effects become important 
is closely correlated with the Kozai resonance. However, the pro- 
cess by which the second most energetic binary migrates inwards 
in the excavated core and sometimes acquires a large inclination 
(near 90°) needs to be investigated further. 

Idealized Kozai cycles may be studied by the three-body reg- 
ularization code which employs the same GR terms. Figure 6 il- 
lustrates the behaviour of the inner eccentricity for a long-lived 
hierarchical triple with initial inclination near 90° using c = 
5 000, 10 000, 15 000 to represent different mass or length scales. 
Simulation mass units are adopted with ao = 1 x 10 -5 , a\ — 
4.5 x 10~ 5 and eccentricities eo = 0.95, e ou t = 0.2. From these 
parameters, e max = 0.9999 while the peak at 0.9980 is reached 
for c = 15 000, with c similar to model W3. This may be suffi- 
cient for GR shrinkage to be achieved on a longer time-scale. In 
this example, the peak eccentricity is maintained over ~ 100 Kozai 
cycles with pericentre advance Aui ~ 0.006 radians per orbit, and 



the slope of the radiation energy loss corresponds to an inspiralling 
time a I a ~ 100. The basic Kozai period is a factor of 5 shorter 
than given by equation (13). Comparison with an earlier analysis 
(Miller & Hamilton 2002, eq. 6) may also be of interest. We ob- 
tain flpN — 0.29 for the dimensionless post-Newtonian Hamilto- 
nian in the realistic example. Although consistent, the expression 
e ma x = 1 — #pn /9 applies to the restricted three-body problem. 

In a second example, ao = 5 x 10 -5 , a\ = 2.3 x 10 -4 , with 
c = 15 000 and inclination 100°. Here the inner semi-major axis 
shrank to ~ 1 x 10~ 5 after some 3 x 10 4 outer periods. During the 
early stage, up to ao — 3 x 10~ 5 , the eccentricity was maintained 
near 0.999 as the result of competition betwen the Kozai cycle and 
relativistic decay. This was eventually superseded by pure decay 
consistent with the analytical expression for e (Peters 1964). 

The case of hyperbolic encounters with the BH has been con- 
sidered in a few simulations only. For this purpose we introduce the 
capture radius (Quinlan & Shapiro 1989), 



= b 



m, mi(m + m i ) 3/ ' 2 1 2/7 



(15) 



where b = 2.68 and «oo is the velocity at infinity in TV-body units. 
With the present parameters in the standard white dwarf model we 



obtain r c 



9 x 10 for a typical excess velocity in the core 



Vac, — 1 and hence capture due to gravitational radiation is not 
likely to occur before disruption. On the other hand, the capture 
process is quite feasible in the point-mass case. 



6 DISCUSSION 

The work presented in this paper should be considered exploratory. 
Although the post-Newtonian formulation itself is well known, a 
practical implementation presents many additional problems. The 
attempt to introduce a scheme based on different time-scales with- 
out sacrificing essential dynamics appears to have been success- 
ful. Likewise, the numerical challenges facing the new implemen- 
tation have been resolved. Thus it has been demonstrated that the 
GRAPE-6 combined with the wheel-spoke regularization code is 
suitable for studying stellar systems with 10 5 members and an 
intermediate massive single BH. A closer analysis of the perfor- 
mance shows that the overheads connected with the central sub- 
system represent less than 10 per cent (for 64-bit operating system 
on GRAPE-6A) which is a small price for an accurate treatment. 
The surprisingly short time-scale for interesting developments also 
means that significant results can be obtained in a few days even 
with the smaller GRAPE-6A. 

One advantage of using a multiple regularization method, as 
opposed to the two-body formulation of nbodyg (Aarseth 2003a), 
is that the complicated derivatives of the GR terms are not required 
(cf. Kupi et al. 2006). On the other hand, it appears that the slight 
defect of introducing a small softening between the low-mass mem- 
bers of the subsystem is justified (and even smaller values may be 
acceptable). So far, the 3PN terms have been added for the short- 
est time-scales, tgr. This was done for illustrative purposes and 
may not be required since the 1PN and 2PN terms act to de-tune 
the Kozai resonance. Contributions from the 3.5PN terms would in 
principle give rise to a recoil velocity; however, this is likely to be 
very small for the present mass ratio. 



2 A scaled softening eo = 0.001 and i? grs 
to a white dwarf radius in model W3. 



10 would correspond 
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So far tidal capture and oblateness effects have not been in- 
cluded (Mardling & Aarseth 2001). General considerations suggest 
that this process would tend to increase the supply of stars into the 
loss-cone. On the technical side, the implementation would be more 
complicated. Thus with two-body regularization it is convenient to 
apply the orbital modifications as an impulse at pericentre while 
in chain regularization this procedure requires more care (Aarseth 
2003a). Likewise, hyperbolic capture of neutron stars or point-mass 
particles by gravitational gravitation (cf. equation (15)) would re- 
quire a similar treatment and may be worth the effort. 

Issues connected with a larger length-scale should also be ad- 
dressed. From equation (7) and the definition of c we have tgr oc 

5/2 

r h in A-body units for the same semi-major axis and eccentric- 
ity. Hence the two-body elements need to be more favourable in 
order for the GR regime to be reached with a comparable effort. 
Increasing the half-mass radii of the white dwarf cluster models 
indicates that relativistic effects may still be appreciable. Further 
explorations of such models are therefore desirable. 

The intriguing question of the relative importance of resonant 
relaxation (Rauch & Tremaine 1996) versus Kozai cycles needs to 
be investigated. Other authors have pointed to the connection (Hop- 
man & Alexander 2006) and, as reported above, there is no doubt 
that the latter process plays a crucial role triggering the early energy 
loss. Typical examples have been identified showing extreme ec- 
centricity evolution for long-lived configurations with semi-major 
axis ratio of about 10. Even so, the actual duration may well be 
too short to be resolved by the present plotting procedure in some 
cases. As for de-tuning of the Kozai cycle, characteristic values 
Alu ~ 0.05 are obtained from equation (14) for the first-order pre- 
cession during critical stages involving white dwarfs and this may 
not be excessive. 

The importance of Kozai cycles in A-body simulations has 
been emphasized before (Iwasawa et al. 2006). In this work, the 
interaction of three massive objects (ttibh = 0.01) in a system of 
N ~ 10 J low-mass particles included the gravitational radiation 
term with fairly small values of c. The formation of hierarchical 
triples involving the three massive bodies frequently led to the ejec- 
tion of one member by the sling-shot mechanism or, alternatively, 
to coalescence of the inner binary. The corresponding eccentricity 
displayed large and small spikes which is characteristic of Kozai 
cycles. As demonstrated above, including the GR precession does 
not suppress the eccentricity spikes. Given the enhanced precession 
rates during the approach to GR conditions, more work is needed 
to clarify this behaviour. In this connection, we note that the pre- 
dicted maximum eccentricity is often reached before triggering the 
final stages leading to coalescence or tidal disruption. 

It is interesting to compare the numerical challenges of study- 
ing one or two massive objects in a stellar system. In the latter 
case, the binary acts to clear the inner region by ejecting stars. The 
depletion of short periods therefore makes for an easier technical 
problem, provided a reliable method is available to deal with the 
massive binary Mikkola & Aarseth 2002, Aarseth 2003b). In either 
case quite large eccentricities are required for GR effects to man- 
ifest themselves. As far as the study of a single BH object is con- 
cerned, the wheel-spoke formulation provides a practical scheme 
for including post-Newtonian terms in a direct A-body code. 

Applications to a star cluster model containing white dwarfs 
also show that significant orbital shrinkage by energy radiation loss 
is possible before disruption. Moreover, high eccentricities by the 
Kozai mechanism can still be achieved. These processes lead to an 
enhanced rate of swallowing by the BH but more extensive simula- 
tions are needed before the growth rate can be determined. Hence 



the present software development provides a powerful tool for ex- 
ploring both idealized and realistic problems of current interest. It 
remains to be seen whether such systems would bear any imprints 
of the swallowing process in the form of supernova events caused 
by white dwarf detonation (Dearborn, Wilson & Mathews 2005). 
Finally, the purpose of the present work is to present a viable new 
method, together with some results illustrating its capability. Hope- 
fully, future work will address issues relating to time-scales and 
model dependence. 
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